Chemo-mechanical instabilities in polarizable active layers 
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We formulate and explore a generic continuum model of a polarizable active layer with neo-Hookean elas- 
ticity and chemo-mechanical interactions. Homogeneous solutions of the model equations exhibit a stationary 
long- wave instability when the medium is activated by expansion, and an oscillatory short-wave instability in the 
case of compressive activation. Both regimes are investigated analytically and numerically. The long-wave in- 
stability initiates a coarsening process, which provides a possible mechanism for the establishment of permanent 
polarization in spherical geometry. 
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Active media are capable to modify their mechanical prop- 
erties under the influence of chemical signaling; in their turn, 
chemical transformations may be affected by local stress, 
strain and polarization. This creates a feedback loop that, on 
the one hand, makes rheology of the medium flexible and re- 
sponsive to the environment and, on the other hand, may cause 
various instabilities, leading to spatial inhomogeneities, oscil- 
lations, and spontaneous motion. The primary example of an 
active biological medium is cy to skeleton, which is the princi- 
pal structural element of eukaryotic cells responsible for their 
integrity, reshaping, proliferation, and motion ||T][2l. 

Experimental studies of living tissues give many examples 
of instabilities involving chemo-mechanical interactions, such 
as spontaneous oscillatory contraction of muscle fibers 131 , 
traveling waves fT], and intermittent localized oscillations in 
vivo that occur, in particular, during such tissue reshap- 
ing processes as dorsal closure (61 |7l and elongation of the 
Drosophila egg shell |8|. Inhomogeneities and patterns of 
non-genetic origin are inconspicuous in vivo but are promi- 
nent in model systems, seen as clusters, swirls and intercon- 
nected bands in motility assays |9|, localized motor activity 
spots in reconstituted actin-myosin mixtures 1 10], and veloc- 
ity spurts in spreading tissues ilTTI . Mechanical impact on 
chemical oscillations (though without feedback reaction) has 
been also detected in inorganic gels 1 12 |. 

Mechanisms of chemo-mechanical interactions in living 
cells and tissues are variegated and largely unknown but un- 
derlying general characteristics can be understood by taking 
as a basis well-established principles of continuous mechan- 
ics that account for the symmetries of the medium, and im- 
posing a simple chemical feedback. The continuum approxi- 
mation is particularly suitable for analytical studies of insta- 
bility thresholds. Symmetry-breaking and oscillatory instabil- 
ities due to chemo-mechanical coupling have been detected in 
this way in the context of calcium transfer dynamics |[T3llT4l . 
and much later with myosin motors acting as chemical vari- 
ables |15|. Further studies included nematic polarization in 
the continuum theory C6l|T3 based on de Gennes' descrip- 
tion of nematic liquid crystals fTSl and adding a term propor- 
tional to chemical potential to account for the biological activ- 
ity. The resulting equations of polar gels are quite formidable 
but their simplified versions have enabled elegant description 



of structures and defects in nematic gels ifTvl, spontaneous 
flow |[T9j[20l, and other phenomena in polar active media. 

In this Letter, we explore an alternative case of a medium 
with vector polarization. The source of polarization of cy- 
toskeleton is directionality of actin filaments characterized by 
a unique direction of treadmilling and motion of myosin mo- 
tors. In acto-myosin stress fibers, nematic order prevails due 
to periodic polarity alternation as, for example, in sarcomeric 
structures. On the other hand, stress fibers in motile cells have 
a graded polarity pattern |21 1 implying vector polarization. 
This kind of polarization is essential for generation of active 
forces that play a primary role in tissue spreading 122^241 . 
Another distinctive feature of our approach is allowing for 
spontaneous emergence of polarization, alongside with me- 
chanical deformation and chemical activity, whereas in earlier 
studies of polar active media, with few exceptions 1251 , only 
the direction, rather than the absolute value of the nematic 
order parameter was evolving dynamically. After formulat- 
ing dynamic equations applicable to a thin flat polarizable de- 
formable layer on a solid substrate, we will determine analyt- 
ically the various instabilities thresholds, and further explore 
numerically the evolution of emerging patterns using the neo- 
Hookean formulation applicable to strongly deformable soft 
matter |26|. 

Dynamics of thin layers dominated by friction is Aris- 
totelian (overdamped), with velocity rather than acceleration 
proportional to the force. Close to equilibrium, the force driv- 
ing the evolution of any dynamic variable ^ can be derived 
from an appropriate energy functional ^{^}, leading to a dy- 
namic equation in the variational form FD^^ = —5J^/5^, 
where = 9t+v- V is the substantial derivative, v is the local 
velocity, and F is an appropriate friction coefficient. We will 
write in this form basic dynamic equations for the deforma- 
tion and polarization fields, and further add non-equilibrium 
terms accounting for the active character of the medium. 

The medium is assumed to be incompressible but it can 
stretch in the horizontal plane x, while modifying the layer 
thickness h. Since the vertical deformation is small compared 
to the deformations in the horizontal plane and the latter, in 
turn, are uniform in the vertical direction, the elastic energy is 
expressed through the two-dimensional (2D) strain tensor. Let 
X, X(x) be the dimensionless position vectors, respectively. 
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in the reference (stress-free) and deformed states. The com- 



ponents of the strain tensor are Ui 



-Sij, where 



5ij is the Kronecker delta; summation over repeated indices is 
presumed here and below. The elastic energy is 



(1) 



where ji is the shear modulus. 

Unlike the standard 2D formulation, the variation of h 
should be taken into account when the dynamic equations are 
derived. In view of the incompressibility, h = hoJ~^, where 
ho is the unperturbed thickness and J is the determinant of the 
2D matrix with the components djXi. We write the dynamic 
equations in the dimensionless form 
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where F = (/i/aa^) f is an active force per unit volume (to be 
specified below), r] = Vqio^/ (2iihQat), a^^ at are the length 
and time scales, and eik is the 2D antisymmetric matrix. 

Planar polarization is characterized by the 2D vector order 
parameter p(X, t). The polarization energy is expressed as an 
integral over the deformed domain: = J i2pd^X, where 
the Lagrangian Cp should be expressed through rotationally 
invariant combinations of p and the 2D gradient operator V 
with respect to the coordinates Xi . Retaining terms with the 
added orders in p and V not exceeding four, we write 
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The second term is irrelevant, unless the coefficient cou- 
ples the order parameter with another field. In an active 
medium, this might be the concentration c of some chemi- 
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cal species. We choose ' as the polarization scale and set 
= —Kq ' c, where c(X, t) is interpreted as a scaled devi- 
ation from a reference concentration cq. Setting also a^s = 0, 
we write the dynamic equation generated by the Lagrangian 
([3]) in the dimensionless form 



V,{J-^c)-a^J-^ (1 + |P|' 



I Pi, 
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where 7 = rpa1/{hzhoat). Generally, one has to include 
the variation of the polarization energy into the deformation 
equations; indeed, a polarized medium may deform to releave 
polarization gradients. We will avoid, however, this compli- 
cation, assuming <C /i. The coupling between polarization 
and deformation will be established instead by the active force 
in Eq. ^ directed along the polarization direction, fi = 2qpi. 

The feedback loop is closed by making concentration of 
the active species strain-dependent. We write the equation of 
the concentration field taking into account that, due to incom- 
pressibility, the concentration is not affected by deformation; 
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FIG. 1: Stability diagram in the parametric plane a, fiq. The oscil- 
latory instability limits are shown for 7 = 5/4 and 77 = 1/3 (solid 
line), ?7 = 5/3 (dashes), and r\ — 10/3 (dots). Inset: the critical 
wavelength k at the oscillatory bifurcation lines parameterized by a. 



the diffusional flux is, however, proportional to the local thick- 
ness. Allowing for the production rate proportional to the lo- 
cal strain and linear decay with the constant identified with the 
inverse time scale l/a^, we write the concentration equation 
in the dimensionless form 



D,c = JV,( V,c) - c + /3( J - 1). 



(5) 



where we have chosen as the length scale the diffusion length 
= VDat, D being the diffusivity. 

The above equation system becomes extremely convoluted 
if applied in a fixed coordinate frame where, besides com- 
puting the substantial derivative that includes the velocity 
V = 9tX, one would need to transform the V operator in 
Eqs. (|4]), (|5j to account for deformations. These complications 
do not affect, however, linear stability analysis, and will be 
avoided at large deformations by implementing a Lagrangean 
numerical procedure. 

We explore now stability of the homogeneous solution p = 
0,c = 0,X = xto infinitesimal perturbations. The linearized 
system determining the bifurcation threshold can be expressed 
through the 2D strain = J — 1 = Vi(X^ — Xi) and splay 
4^ = ^iPi by taking the divergence of linearized Eqs. ([2]), (|4]): 
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A nontrivial solution of Eqs. ([6]) - ([5]) is sought for as a wave- 
form exp(At + i k • x) with a wave vector k. The eigenvalues 
A determining stability of the uniform state c = = (j) = 
are computed by setting to zero the determinant of the Fourier 
transformed linearized system. The bifurcation thresholds are 
computed in a standard way by requiring sup|p,| (3? A) = 0. 

In the case of extentional activation (/3 > 0), one can show 
that the branch of real eigenvalues first reaches zero in the 



FIG. 2: Snapshots from a simulation in the oscillatory unstable para- 
metric region (/3q = —210, {I3q)c ~ —196). The time difference 
between the two images equals half the oscillation period. The sim- 
ulation domain is a double-periodic 10 a/5 x 10 a/5 square. Lighter 
shades correspond to higher concentration levels; scaled arrows indi- 
cate the polarization p. 

long- wave mode k = at /3g = a^. In the case of com- 
pressive activation (/3 < 0), a wave instability at finite |k| and 
7^ is observed. The analytic formulas for the critical 
wavenumber kc and the critical value /3c are too cumbersome 
to be of practical use and will not be listed here. The min- 
imal gap width between the two instability limits equals to 
7 + 7^. The bifurcation diagram in the parametric plane a, Pq 
is shown in Fig.[T] 

We explore dynamics beyond the bifurcation thresholds for 
both types of activation by direct numerical simulation, start- 
ing from randomly perturbed homogeneous initial conditions. 
We discretize Eqs. ([2]), (|4]), and ^ using a flexible grid with 
nodes following the physical displacement of the medium. In 
the spirit of mesh-free computations |27|, the spatial deriva- 
tives are calculated based on local approximation of the fields 
by low-order polynomials. The resulting set of ordinary dif- 
ferential equations is then solved using an embedded adaptive 
Runge-Kutta-Fehlberg scheme of the order 4 and 5 |28|. The 
parameters are set to q = a/5/100, a = a/5, = = 
5/3, 7 = 5/4. 

In the oscillatory regime, we observe a wave pattern with 
the active force (directed from high to low c regions) follow- 
ing the rotating concentration gradient between the domains 
oscillating in antiphase, as seen in Fig. |2] The directions of 
the elastic and active force are, however, disaligned due to de- 
lay, so that the angle between the two is widely distributed 
with a shallow maximum about 7r/2 (Fig.jsJ. 

In the case of extensional activation, the behavior is qual- 
itatively different. Initially, a short-wave pattern correspond- 
ing to a fastest-growing mode emerges from the initial noise. 
Then, a slow coarsening process ensues, during which larger 
domains grow while smaller ones either merge or vanish, so 
that the dominant wavelength grows, as seen in Fig. |4] The 
wavelength corresponding to the maximum of the Fourier 
spectrum is plotted against time in Fig. [5] Due to the finite 
size of the simulation domain, the plot for each individual run 
has a number of discontinuous steps but it approaches a mono- 
tonically growing curve when averaged over 30 runs. 



FIG. 3: Histogram of the angle cp between the directions of the 
elastic and active force. 




FIG. 4: Domain coarsening: snapshots at t = 4 x 10^ (a), 2 x 
10^ (b), 3 X 10^ (c), 5 X 10^ (d), simulation on the double-periodic 
100^5 X 100^5 square domain for /3q = 6, {/3q)c = 5. 
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FIG. 5: Time evolution of the wavelength of the dominant Fourier 
mode during coarsening (the average of 30 simulations on the 
double-periodic 100^5 x 100^5 square domain for I3q — 6). 
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FIG. 6: Cross section through the dashed Hne in Fig. |4] (d). The 
polarization component along the line p • and deformation J — 1 
are scaled to show all three fields in the same plot. 
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FIG. 7: Time evolution of the concentration c at different polar angles 
^ on a sphere of radius 10 a/5. 



Unlike the oscillatory case, the active and elastic forces are 
almost perfectly anti-aligned: the mutual angle deviates from 
TT by less than 3.5° with 98% probability. As the active force 
is directed against the concentration gradient, the low-c re- 
gions are squeezed into narrow stripes separating high-c do- 
mains. In turn, the chemical production is inhibited within 
these stripes due to strong compression, as the deformation 
acts as the source of c in Eq. ([5]). The resulting correlation 
between J and c is seen in Fig.[6]that shows the cross-section 
profiles at a late stage of coarsening. 

The coarsening process provides a natural mechanism for 
polarity establishment in cortical tissues, as is observed, for 
example, in the C. elegans embryo 1291 . When Eqs. ([2]), 
(|4]), and ([5]) are solved on a spherical surface, the coarsen- 
ing consolidates two domains, whereafter the smaller one is 
compressed into a low-c spot. This is illustrated in Fig. [7] 
where the evolution of c is plotted for five different polar an- 
gles 6, where 6> = is defined by the location of the high- 
concentration pole. As the polarization is aligned with the 
concentration gradient, the active force is directed from one 
pole to the other. The established polarity axis is permanent, 
unlike Ref. |[30ll where advection serves as a mere trigger re- 
quiring additional means for preserving the polarization. 
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